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We study the phase diagram of a three-component Fermi gas with weak attractive interactions, 
which shows three superfluid and one normal phases. At weak symmetry breaking between the 
components the existence of domain walls interpolating between two superfluids introduces a new 
length scale much larger than the coherence length of each superfluid. This, in particular, limits 
the applicability of the local density approximation in the trapped case, which we also discuss. In 
the same regime the system hosts soft collective modes with a mass much smaller than the energy 
gaps of individual superfluids. We derive their dispersion relations at zero and finite temperatures 
and demonstrate that their presence leads to a significant enhancement of fluctuations near the 
superfluid-normal transitions. 
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I. INTRODUCTION 

In Landau's approach to phase transitions, conven- 
tional superfluidity and superconductivity are charac- 
terized by a single complex order parameter. However, 
in certain instances a proper description of a superfluid 
within this approach requires the introduction of an or- 
der parameter with several complex components. For 
example, in the case of superfluidity in '^He [1] the spin- 
triplet, p-wave pairing is described by nine coefficients, 
and three different superffuid phases are experimentally 
realized. Other examples include unconventional super- 
conductivity [2] in heavy fermion coinpounds and color 
superconductivity in nuclear matter [3|, where different 
phases could be realized depending on the chemical po- 
tentials - the two-flavor color superconductor and the 
color-flavor-locked phase. Many aspects of multicompo- 
nent superfluidity in these systems are understood only 
on a phenomenological level due to their intrinsic com- 
plexity. Atomic Fermi gases, on the other hand, provide 
a unique avenue to explore these phenomena in a highly 
controllable way, thanks to the tunability of the interac- 
tions between atoms. Multicomponent superfluidity in 
atomic fermions could be realized, for example, by trap- 
ping and cooliiig multiple hyperfine states of the same 
atomic species 0, Q or of different species Q . 

Here we present a study of the phase diagram of a 
three-component Fermi gas with weak attractive interac- 
tions. In particular, we consider the situation in which 
the "color" symmetry between the components is broken 
due to differences in the interaction strengths or chemical 
potentials, while the masses are the same. This situation 
is relevant to possible experiments involving three hy- 
perfine states of ^Li atoms. We first develop a Ginzburg- 
Landau expansion for this system and use it to confirm 
the previous results 0, H, that there are four possible 
phases: the normal state and three superfluid states Si, 
S2, and S3. In each of the superfluid states two out of 
three components are paired, while the third one is in a 



normal state. First order phase transitions between dif- 
ferent superfluid states can be driven by varying interac- 
tion strengths, chemical potentials, particle densities or 
temperature. We construct the phase diagram in grand 
canonical and canonical ensembles at finite and zero tem- 
peratures, shown in Figs. [T]- [31 and [HI see also Refs. [8, ^. 

The canonical phase diagram at fixed temperature 
(Fig. [2) contains regions where the homogeneous state 
is unstable. When the particle densities are within these 
regions the two superfluid states phase separate, as ex- 
pected for the first-order phase transition; see, e.g., p^ . 
We therefore explore the properties of domain walls be- 
tween different superfluids; see Fig. |4l In particular, we 
explicitly determine the shape and the thickness i of the 
domain walls in various regimes. The case of weak sym- 
metry breaking between two components (say, 1 and 2) 
- i.e., when the couplings of 1 and 2 with 3 and chemical 
potentials /ii and ^2 are close - is especially interesting. 
At full symmetry between 1 and 2, £ — 00. This is natural 
as in this case the thermodynamic potential r2(Ai, A2), 
where Ai and A2 are the order parameters for superfluid 
states Si and S2, respectively, is invariant with respect to 
rotations in the A1-A2 space. This implies that the two 
minima of the potential (Ai, A2) = {(A?, 0), (0, A§)} de- 
scribing superfluids Si and S2 can be connected by con- 
tinuous lines of minima. Then, Ai can be continuously 
deformed into A2 at no energy cost when moving from 
one point in space to another. At weak symmetry break- 
ing, as we demonstrate below, the thickness of the do- 
main wall £ is parametrically larger than the coherence 
lengths ^1 and ^2 of superfluids Si and S2. 

For a trapped three-component gas the local density 
approximation (LDA) predicts sharp boundaries between 
superfluid states Si and S2 [IH. In reality, there is a do- 
main wall of length £ between Si and S2 where the two 
superfluids coexist. Therefore, the characteristic length 
scale over which the boundaries predicted by the LDA 
are smeared is £, rather than ^1 or ^2- Moreover, if the 
radius of the trap, R, is comparable to £, the two su- 
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perfluids coexist throughout the trap, so that the cases 
£ and R ~ d are quahtatively different. In particu- 
lar, this means that the LDA breaks down in the entire 
trap when R ^ i. For typical experimental parameters, 
the condition R ^ £ translates into Nt ^ 10^ (see be- 
low), where Nt is the total number of fermions of all three 
species. Moreover, as we will see, the deviations from the 
LDA are significant already for Nt as large as 10^. 

Another consequence of the weak symmetry breaking 
discussed above is the presence of soft collective modes 
in multicomponent Fermi gases 0, Suppose, for 

example, the system is in the superfluid state Si. By 
considering the thermodynamic potential f7(Ai,A2) as 
above in the case of the domain wall, we expect fluctu- 
ations i5A2(r,t) towards superfluid S2 to be massless in 
the symmetric case. Below we derive the mass and dis- 
persion relations of the corresponding collective modes in 
the general asymmetric case at T = and at finite tem- 
peratures. We show that at weak symmetry breaking the 
mass can be much smaller than the BCS gap in the su- 
perfluid Si in the absence of the third fermionic species. 
At finite temperature these fluctuations result, in partic- 
ular, in an enhancement of the Ginzburg-Levanyuk num- 
ber Gi by a large factor. In other words, the window 
of temperatures around the critical temperature for the 
normal - superfluid Si transition where fluctuations dom- 
inate becomes much larger in the presence of the third 
component. 

Let us comment on the experimental realization of su- 
perfluidity in a three-component Fermi gas. Achieving 
a stable gas in this case appears more challenging than 
in a two-component one due to the enhanced role of the 
three-body scattering. In the two-component Fermi gas 
three-body recombination is suppressed thanks to the 
Pauli exclusion principle fl^ and the system is stable 
over tens of seconds. In the three-component case there 
is no such suppression and the decay time is of the or- 
der of milliseconds 0. Various proposals are being put 
forward in order to increase the lifetime of the system, 
such as, e.g., the stabilization by an optical lattice [T^ 
similar to that for bosonic atoms We note that the 
results presented in this paper are obtained in the weak- 
coupling regime, which is expected to be insensitive to 
the stabilization technique. For example, a lattice added 
to the trapping potential affects the single-fermion spec- 
trum only. This is irrelevant at weak coupling since the 
superfluid energy scales are assumed to be much smaller 
than the fermionic bandwidth. The single-particle bands 
contribute only through the density of states at the Fermi 
energy irrespective of the details of the spectrum. 

The paper is organized as follows. In the next sec- 
tion we give a brief overview of the mean-field approach 
and introduce our notation. In Sec. IIIII we present the 
Ginzburg-Landau expansion of the thermodynamic po- 
tential and discuss the phase diagram at finite temper- 
atures. We study the domain walls in Sec. |TVl and in 
Sec.|V]we describe the zero-temperature phase diagram. 
Section IVII is devoted to the collective modes. Finally, 



we summarize our results in Sec. I VIII 



II. THERMODYNAMIC POTENTIAL 

In this section we outline the derivation of the thermo- 
dynamic potential, from which the phase diagram and all 
thermodynamic quantities can be obtained. We will not 
go into details, as the derivation is a well-known proce- 
dure [iBl- Our starting point is the following Hamilto- 



(1) 



Here Hq — p^/(2m) is the single-particle Hamiltonian 
(we assume that all the particles have the same mass). 
As discussed in the Introduction, an optical lattice would 
modify the single-particle Hamiltonian. Its effect can 
be taken into account by introducing an effective mass 
"T'cff 7^ 'm-, which in the weak-coupling regime results only 
in a renormalization of the density of states introduced 
below in Eq. ([7]). The pairwise interaction part is 
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(-^jeyfeV'i) {il^k'£t]'k''4'j') , (2) 



where Sijk is the totally antisymmetric tensor and 
{i, j, fc} = {1, 2, 3}. By the Hubbard- Stratonovich trans- 
formation, we introduce the pairing field A(t, r) = 
(Ai(t, r), A2(t, r), A3(r, r)) and after integrating out 
the particle fields we obtain the following effective 
action for A(r, r): 
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drdrx 
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(3) 



where = diag (5^^^) and G^^ is the particles' in- 
verse Green's function, which is a 6 x 6 matrix in Nambu- 
Gorkov space with the structure 



G- 



(4) 

Here /ii are the chemical potentials for the different 
species. 

In the mean-field approximation the thermodynamic 
potential is obtained by evaluating the effective action ([3]) 
for a T-independent pairing field A(r). This is expected 
to be an excellent approximation for the description of 
a weakly coupled fermionic superfluid at temperatures 
not extremely close to the transition temperature p^ . 
First, let us consider the case of a uniform order param- 
eter. Performing a Fourier transform from real space- 
imaginary time to the momentum-Matsubara-frequency 
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space in Eq. ([3]) we derive 



(27r)3 



p 

+ ^ lAiplAjf [ (w„ + i^i) (ujn - i^j) + c.c. 



(5) 



where aj„ = 2TTT{n + 1/2), = p^/(2to) — fit, the 
sum over P denotes the sum over cychc permutations 
of {i,j,k} ~ {1,2,3}, and "c.c." is the complex conju- 
gate. For vanishing order parameter A = 0, we obtain 
the sum of the thermodynamic potentials for three per- 
fect g expected. Also, for an order parameter 
with only one nonvanishing component A.^ ^ 0, Eq. ([5]) 
reduces to the sum of the potentials of a normal gas 
and a two-component Fermi superfluid. Let us denote 
the corresponding zero-temperature order parameter of 
the two-component superfluid in the absence of the third 
fermionic species as A". We note that Eq. ([5]) is ultra- 
violet divergent, and a regularization procedure (e.g., a 
hard cutoff as for superconductors 'To] or a T-matrix ap- 
proach [l3|) should be implemented. Then all physical 



III. GINZBURG-LANDAU EXPANSION 

Here we perform a Ginzburg-Landau expansion for the 
thermodynamic potential and use it to obtain the finite- 
temperature phase diagram of the system in the /ii-/i2 
plane; see Fig.[TJ We determine the superfluid-superfluid 
and superfluid-normal transition lines and the metasta- 
bility regions in both grand-canonical (Fig. [T]) and canon- 
ical (Fig. [2) ensembles. In the latter case there is a re- 
gion of the phase diagram where a homogeneous state 
is unstable and a phase separation between two types of 
superfluid takes place. We identify this region as well as 
the corresponding supercooling lines; see Fig. O 

According to Landau's phenomenological approach 
[Toj . the thermodynamic potential near a second-order 
phase transition can be expanded in powers of the or- 
der parameter. If only even powers are present and the 
coefficient of the fourth-order term is positive, the van- 
ishing of the coefhcient of the quadratic term determines 
the second-order transition point. When the fourth-order 
term also changes sign, the transition becomes first order, 
and higher-order terms should be included in the power 
series. 

As shown by Gorkov [3, this phenomenological the- 
ory can be derived by expanding the microscopic theory 
in |A|/27rr around A = 0. Using Eq. © [or Eq. ® 
for the uniform part], we obtain to the fourth order in 
components of A = (Ai, A2, 0) 



do in what follows. 

The (meta) stable states are given by the (local) min- 
ima of Q. This condition determines the mean- field phase 
diagram. We will show below that the possible phases fall 
into two classes - normal state or a two component su- 
perfluid plus a normal gas - in agreement with the results 
of [i, Q . The two superfluid components can be any two 
of the three atomic species; i.e., there are three possi- 
ble superfluid states, which we denote as Si, S2, and S3 
when the paired species are 2 and 3, 1 and 3, and 1 and 
2, respectively. For simplicity, unless otherwise specified, 
we assume from now on 173 = 0, so that A3 = and 
only Ai and A2 components of the order parameter can 
be nonzero. This can be a good approximation in the 
case of three hyperfine states of ^Li, where two out of 
three Feshbach resonances mediating the attractive in- 
teractions between the states are close in magnetic fields 
[J . The third resonance is at a lower field and can be ne- 
glected on the BCS side of the crossover. The inclusion of 
the case 173 ^ in our formalism is straightforward. We 
briefly comment on this case in Sec. IIII CI and show the 
corresponding phase diagram in Fig. [31 For concreteness, 
we take > \g2\ and introduce the notation: 

/ii = M3 - Ai2 , /i2 = - Ml : (6) 
for the differences in chemical potentials. 



quantities can be expressed in terms of the A°'s, as we fl — Hn ='^E/ + ~ 



|A,|4 + i|:|VA,p 



i./3i2|AinA2p, 



(7) 



where i7Ar is the normal-state thermodynamic potential 
of the ideal gas, ly is the density of states at the Fermi 
energy, and the coefficients a^. Pi, and /?i2 are 



a, ^In— +Re* (-+[J^) - ^ (- 
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and 
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hi - h2 AttT 



-Im 



2 4ttT / ' 



\2 AttTJ 



(8) 



(9) 



2 47rr 



(10) 



Here '^{x) is the digamma function and is the critical 
temperature of the superfluid Si at zero chemical poten- 
tial difference and in the absence of the third fermionic 
species - i.e., for hi = and gj^i = 0. According to the 
standard BCS theory for two species, T^^ is related to 
the corresponding zero-temperature order parameter A? 
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as Tc- = c'^^A^/tt, where 7^ is Euler's constant. Note 
that due to our choice \gi\ > \g2\ for the coupling con- 
stants, Tcj^ > Tc2- We will comment below on the physi- 
cal meaning of the temperatures Tc- in the three species 
case. 

Expressions ([8]), and (fTO|l for the coefficients in 
the Ginzburg-Landau expansion ([7]) were derived in the 
weak-coupling limit, which enabled us to approximate 
the density of states with a constant v. However, the 
structure of the potential ([7]) is dictated by symmetry and 
must remain the same at any coupling. Indeed, since the 
particle number is conserved separately for each species, 
the potential must be independent of the phases of the 
complex components of the order parameter. Therefore, 
the only allowed terms in the expansion to the fourth 
order are |A,|2, |A,|4, and lAipjAap. 

The terms in curly brackets in Eq. ([7]) give the ther- 
modynamic potential Oi(Ai, hi, T) of the two component 
superfluid Si in the absence of the third species, while the 
|Aip|A2p term represents the interaction between the 
two superfluids. The same expression for fli{Ai, hi,T) 
was previously obtained (Toj in a study of the nonuni- 
form superconducting Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) state In thin superconducting films in 

a parallel magnetic field the thermodynamic potential 
ili{Ai, hi, T) describes the effect of the Zeeman splitting. 
In this case, hi has a meaning of the Zeeman magnetic 
field and A^ is the superconducting order parameter. Let 
us briefly summarize the phases described by fli in the 
hi-T plane [2l[ before we proceed to the phase diagram 
for three species. For T > Tc- the quadratic coefficient 
is positive, > 0, and the two-component Fermi gas 
is in the normal state for any value of hi. At temper- 
atures Tjj.; < T < Tc- a. second-order transition to the 
superfluid state Si takes place when ai(hi,T) — 0. For 
temperatures lower than the tricritical temperature. 



OMTc, 



1,2, 



(11) 



the quartic coefficient is negative, f3i < 0, whenever 
ai — > and the normal-superfluid Si transition is flrst 
order. The tricritical temperature and the correspond- 
ing tricritical chemical potential are determined from the 
condition a^{hl^,T^^.) = /3,(/ij^i, T^,.) = 0. This picture 
can also be obtained in the BCS limit from the phase dia- 
gram for polarized Fermi gases in the BCS-BEC crossover 

Mm. 

Now we turn to the analysis of general properties 
of the full thermodynamic potential for three species. 
The Ginzburg-Landau expansion ([7]) is a good approx- 
imation only when the polynomial /3i|Ai|^ -I- /32|A2|^ + 
2/3i2|Aip|A2p is positively defined. Otherwise, 
—00 as |A| —^ 00 along a certain direction in the A1-A2 
plane. Using Eqs. (O and pl}|) . one can show that this 
condition reduces to 



11 „ / 1 
Re*" - 

A(2TrT)^ \2 



AttT 



>0, 



1,2. 



These inequalities are equivalent to Ihil/inT < 0.304. 
Since Eq. ([7]) was obtained by expanding in \A\/2ttT, we 
also should have \A\/2ttT ^ 1. Thus, the conditions of 
applicability of the expression ([7]) for the thermodynamic 
potential are 



AttT 



< 0.304, 



2ttT 



« 1, 



1,2. 



(13) 



It follows from the discussion in the previous para- 
graph that these inequalities hold only in the "high- 
temperature" regime T > T^^^ > T^^^. Otherwise, /3i < 
when Ai/2ttT ^ 1 and the first condition in Eq. 
is violated. For the remainder of this section we restrict 
ourselves to this range of temperatures. Then, one can 
show using Eqs. ^ and Q that the condition (|12p al- 
ways holds whenever any of the quadratic coefficients ai 
is sufficiently small. 

Let us discuss the possible phases of the three species 
system in the /ii-ft.2 plane as a function of temperature 
going from higher to lower temperatures. For T > Tc^ > 
Tc2 we see from Eq. ^ that both quadratic coefficients 
in Eq. ([7]) are positive, i.e.. 



, T fl h, 

a.=ln-+Re* 



>0, 



(14) 



(12) 



for any h-i and i = 1,2. In this case, the only stable state 
is A = - i.e., the normal state. As the temperature is 
lowered, ai first vanishes at T = and hi = 0, while a2 
remains positive. Therefore, a phase transition from the 
normal to the superfluid state 5*1 occurs and Tci is the 
actual critical temperature for this transition. Generally, 
for <T < Tc^ , we have a2 > 0, while ai changes sign 
at a temperature-dependent critical chemical potential 
h^T) determined by the equation ai {h1{T),T) = 0. At 
hi < h\{T) the superfluid state 5*1 is the stable one, while 
at larger hi the system turns normal. 

The case T < Tc^ is more complicated. Now the con- 
ditions hold for both components only when both 
\hi\ and |ft,2| are sufficiently large. In the /ii-/i2 plane 
Eq. (fT4|l determines four normal-state regions; see Fig.[T] 
A second-order phase transition from the normal to a su- 
perfluid state takes place when one of the coefficients 
ai changes sign. For example, starting from the normal 
state, keeping /i2 fixed, and changing hi, we get a transi- 
tion between the normal state and the superfiuid Si, as 
ai becomes negative while a2 is still positive. This ar- 
gument, however, cannot predict the state of the system 
in the central region of the /ii-ft.2 plane where both a^'s 
are negative. We will explore this region in more detail 
in the following subsection. 



A. Phase diagram in the vicinity of critical 
temperatures 

As discussed above, the normal state is the stable phase 
in four sectors of the phase diagram, corresponding to the 
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four corners in Fig. [T] Here we show that in the central 
region two different cases are possible: (i) the thermody- 
namic potential has only one minimum, which coincides 
with the superfluid state Si for one of the two possible 
condensates A^; (ii) Q has two local minima, such that 
one condensate is the stable state and the other one is 
a metastable one. In the latter situation, a first-order 
phase transition separates the two superfluid states, as 
identified by the dashed lines in Fig. [TJ The two minima 
are degenerate along these lines in the /ii-/i2 plane. The 
gray areas around the lines shown in Fig. [T] enclose the 
regions where two local minima are present. 

In this subsection we obtain the phase diagram for the 
case when the two coupling constants gi and g2 are suffi- 
ciently close in magnitude. We also take the temperature 
to be near the critical temperatures and Tcj, i.e., 



(15) 



The first inequality in Eq. (fTS]) holds since T^^ is the 
critical temperature for the two component superfluid 
with coupling gi [see the text below Eq. (fTO|) ] and the 
couplings are close. As we will see below, in this case the 
condition |A|/27rT ^ I for the validity of the Ginzburg- 
Landau expansion holds. Then, it follows from Eq. ([T3|) 
that expression ([7]) for the thermodynamic potential can 
be used not just near the phase transition lines, but for all 
hi and /12 such that |/ii|/47rT < 0.304. Nevertheless, the 
conclusions we draw regarding the phase diagram have 
general validity at sufficiently high temperatures T > 
T^i; ~ O.SGTc,; see the text below Eq. ^ and at the 
end of this subsection. 

Let us first consider a homogenous system; i.e., the 
gradient terms in Eq. ([7]) vanish: 



(16) 



To find the stationary points of 51, we pass to a polar 
coordinate representation 



A cos 61, IA2 



A sin 6* . 



(17) 



Differentiating Eq. (jl6[) with respect to the angular vari- 
able 0, we find 



A^ cos 6 sin ( 



Q!2 — ai 



+ A2 ((/3i2 - /3l) C0S2 -f (/?2 - /3l2) sh 



(18) 



where a^, /3i, and /3i2 are defined by Eqs. ([HI), ®, and 
PU)) . respectively. Equation always admits the three 
solutions A = 0, = 0, and 9 = -k/2. These are, respec- 
tively, the normal state, the condensate Ai, and the con- 
densate A2. To determine the value of the nonvanishing 
order parameter component, we also need to equate to 
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FIG. 1: Finite-temperature (T = O.SSTca) phase diagram for 
a three-component Fermi gas in the /11-/12 plane of chemical 
potential differences, Eq. The two nonvanishing pairwise 
couplings (71,2 are such that Tc^/Tc^ ~ 1.04, where are 
defined below Eq. (|lUp . Note that the normal state (N) is 
stable at large hi while the superfluid states (Si) at lower 
ones. The stronger interaction {\gi\ > |(/2|) determines the 
superfluid state (Si) realized at hi — h2 — 0. The horizontal 
(vertical) segments denote second-order N-Si (N-S2) transi- 
tions. The dashed curves mark the first-order S1-S2 transi- 
tions; see Eq. (|23p . The shaded areas limited by the dotted 
curves [Eq. (I22|l ] are the metastability regions. 



zero the derivative of the thermodynamic potential ([T 
with respect to A: 



= A 



ai cos^ 6 + a2 sin^ ( 



A^ cos^ 9 + 132 sin^ 9 + 2/3i2 cos^ 9 sin^ ■ 



(19) 



We obtain 



for 6* = and 



A? 



= 0, 



-02/ 1^2 



(20) 



(21) 



for 9 = 7r/2. We see that A/(27rT) is indeed small near 
the second-order N-Si phase transition, since — > 
at the transition. This also implies that in the super- 
fluid state Si the chemical potential difference is such 
that \h^\ <C 47rT, because the condition \T - < Tc,, 
Eq. (fTS)) . makes the first term in the definition ([8]) of ai 
small. 

Having found the stationary points, we must check 
their stability. For a minimum, the second derivative 
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must be positive. The second derivative of the thermo- 
dynamic potential at stationary points (PU)) and (PT|) 
vanishes when 



aiPj - ajf3i2 = {i ^ j). 



(22) 



These equations define the stability lines enclosing the 
regions with two minima (gray areas in Fig. [1]) . Outside 
these regions, there is only one minimum, while the other 
stationary point is a saddle. Note that the fourth solution 
to Eq. (fT8|) can be obtained by equating the terms in 
square brackets to zero. This solution is present only 
when the thermodynamic potential has two minima and 
corresponds to the saddle point between them. 

Finally, the minima (|^D|) and (|2ip are degenerate 
(the dashed hues in Fig. [T]) when Q{—ai/ Pi,0) = 
il{0, — Q!2//32), which yields 



Pi 02 



(23) 



This condition can be satisfied only close to both second- 
order N-Si transitions, so that ^ AttT for both i — 1 
and 1 = 2; see the discussion after Eq. ([2T|l . Substituting 
Eqs. ^ and ^ into Eq. we obtain to leading order 
in \h^\/4TTT 



1 



hi -hi 
(47rT)2 



= ln^ 



(24) 



This equation shows that chemical potential differences, 
temperature, and the asymmetry in the interaction 
strengths determine the lines of the first-order phase 
transitions between different superfluid states. At fixed 
T, Eq. (|24p defines transition lines hi(h2) in the hi-h2 
plane; see the dashed curves in Fig. [T] 

In the presence of a trapping potential V{r), we can 
combine our phase diagram of Fig. [T] with the so-called 
LDA T?! to predict the formation of different superfluid 
shells in the trap. The LDA assumes position-dependent 
chemical potentials 



^l,=^l'!-Vir), 1 = 1,2,3. 



(25) 



The differences hi = /Xg— /ij and /12 = ^3—^1 [see Eq. ([5])] 
remain constant throughout the trap and identify a point 
h = {hi, h2) on the phase diagram (Fig. [1]). The temper- 
ature = Tc-{^i) depends on the chemical potential 
Hi as in the standard BCS theory; see the text below 
Eq. pn)) . As /ii decreases from the center to the edge of 
the trap, Tc^{fJ,i) also decreases. On the other hand, the 
positions of the lines in the phase diagram in Fig. [T] are 
determined by the values of Tc-{^J,i); see, e.g., Eq. ([M]). 
Therefore, the "local" phase diagram - i.e., the phase 
diagram of the homogenous system that corresponds to 
the values of chemical potentials at a particular point 
r in the trap - changes, and as we move from its cen- 
ter towards the edge, the regions where the superfluids 
are stable become smaller due to the decrease in Tc-{fJ.i). 



The actual values of fJ^ and consequently hi must be de- 
termined self-consistently by fixing particle numbers for 
species 1,2, and 3. Depending on the position of the re- 
sulting point h in the local phase diagram at the trap 
center, different configurations are possible. For exam- 
ple, if h is in the Si region at the center, the evolution 
of the local phase diagram with the position r can bring 
this point into the N region or make it pass through the 
S2 region first. These two possibilities correspond to a 
central superfluid Si core surrounded by a normal shell or 
a superfluid Si core followed by an S2 shell and a normal 
shell farther out, respectively. If h is in the S2 region at 
the trap center, on the other hand, we obtain a superfluid 
S2 core surrounded by a normal shell. Alternatively, for 
low particle number the normal-state atoms of the non- 
condensed species could form a normal core overlapping 
with the superfluid one. This qualitative picture is in 
agreement with numerical results of [11]. However, as we 
will discuss at the end of Sec. IIVI the LDA has rather 
limited applicability in the presence of an S1-S2 bound- 
ary. 

Let us summarize our observations so far in this sec- 
tion about the possible phases and phase transitions in 
the homogeneous case. We saw that for T > Tc^ > Tc^ 
the system is in the normal state N for any chemical po- 
tentials differences hi = ^3 — ^2 and /12 = Ms ^ Mi (recall 
that we set the coupling constant between species 1 
and 2 to zero, while \gi\ > |(?2|)- A second-order phase 
transition to the superfluid state Si where species 2 and 
3 condense first happens at /ii = and T = Tc^ at ar- 
bitrary ft,2- For < T < the only possible states 
are the normal state and superfluid Si. At lower tem- 
peratures T < three states can exist as shown in 
Fig. [1] A second-order transition from the normal state 
to superfluid S2 first takes place at T = T^^, ft-2 = 0, and 
sufhciently large \hi\ (so that S2 wins over Si); see Fig.[T] 
The S1-S2 transition is always first order, while the N-S2 
and N-Si are both second order provided that the tem- 
perature is above the tricritical temperatures (jlip. i.e.. 



T > ma.x{Tl^i,Tli} = Tu 



(26) 



For T < Ttii at least one of the transitions N-Si or N-S2 
becomes first order and the Ginzburg-Landau expansion 
([7]) breaks down; see the discussion below Eq. (|13p. 



B. Phase separation 

In the previous subsections we analyzed the phase dia- 
gram in the grand-canonical ensemble. Here we consider 
the canonical ensemble; i.e., we fix the densities of the 
three fermionic species. 

The corresponding chemical potentials arc found by 
solving the equations 



Oil 
dfii 



(27) 
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If the differences between the densities are large, the 
chemical potential differences are also large and the gas 
is in the normal state. Let us assume that the densities 
deviate little from an average density no, 



Ui = no + Srii 



(28) 



Then, density deviations can be written as the sum of a 
noninteracting term and a correction due to the presence 
of the superfluid: 



dsn 



(29) 



where — fii — fiQ, SQ = Q — flj^, and /iq is the chemical 
potential for a noninteracting gas with density hq. In 
Eq. (|29|) we neglected finite-temperature corrections to 
the noninteracting contribution [l^] • 

For example, if the system is in the superfluid state Si, 
we find using Eqs. H ^ . P^ . and PU)) : 



where 



hi — 713 — 12 = vhi — 2 



2pi 



dhi 



(30) 



(31) 



Using similar equations for the homogenous superfluid 
S2, we obtain the phase diagram presented in Fig. [2]by 
mapping the lines in the phase diagram in the /11-/12 
space of Fig. [T] onto the corresponding lines in the ni- 
n2 space of density differences. In particular, we note 
that each first-order phase transition line in the upper 
and lower half planes of Fig. [T] (dashed lines) maps into 
two lines. Indeed, according to Eqs. ((3T|) and (p9| . Sil 
and therefore hi and n2 are different on the two sides of 
the transition. This means that, as usual in the case of 
first-order phase transitions, there is a region in the phase 
diagram where no homogeneous state is stable and phase 
separation must occur. Between this region and the sta- 
ble homogeneous superfluid states, there are supercool- 
ing regions (gray areas) where the homogeneous states 
are metastable towards phase separation. The limits of 
these regions are found by mapping the corresponding 
limiting metastability lines in the grand-canonical phase 
diagram (dotted curves in Fig. [T]). 

At the end of the previous subsection we argued that, 
within the LDA, the two superfluid states can coexist in 
a trap. In this section we have shown that the transition 
between the two superfluids is necessarily accompanied 
by a jump in the density. In this sense, it is similar 
to the low-temperature transition between the superfluid 
and normal states in the polarized two-component gas 
[l7| . In this system, the density jump signals a potential 
breakdown of the LDA on the length scale of the coher- 
ence length. There is also evidence that surface tension 
effects should be taken into account to explain the shape 
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FIG. 2: Finite-temperature (T = O.SSTca) phase diagram for 
a three-component Fermi gas in the plane of particle density 
differences ni-n2, Eq. pop . As in Fig. [T] the two nonvan- 
ishing pairwise couplings are chosen so that Tci/Tcg ~ 1.04. 
The horizontal (vertical) segments denote the second-order 
phase transition N-Si (N-S2) between normal (N) and su- 
perfluid state Si (S2). Dashed curves represent the limits 
of stability for the homogenous superfluids and enclose the 
phase-separated (PS) states. The shaded areas are the su- 
percooling regions where a homogeneous superfluid state is 
metastable toward phase separation. 



of the superfluid core in elongated traps [25|. In the 
present case of the S1-S2 transition there are two com- 
peting length scales (the two coherence lengths) that can 
affect the properties of the interface, which we study in 
the next section. 



C. Phase diagram in the case when all three 
couplings are nonzero 

Here we briefly discuss the case when the coupling con- 
stant gs between species 1 and 2 is also nonzero. Let 
1 53 1 < Iff2|- Now, in addition to Tc^ ^ there is the third 
temperature scale Tc^ . Similarly to ^ , it is deflned 
as the critical temperature of the superfluid with com- 
ponents 1 and 2 in the absence of 3. Further, additional 
terms containing IA3P appear in the thermodynamic po- 



tential p6)) . The coefficient of IA3P is defined by 
Eq. ([5]) with = |J■2—^J■l = h2 — hi. For T > Tc^, we have 
as > and the phase diagram in Fig. [T] is unchanged. For 
Ttii < T < Tc^ , new N-S3 second-order phase transitions 
are possible as well as first-order transitions S3-S1 and 
S3-S2. A phase diagram with these transitions is shown 
in Fig. [3] Note that, since /13 = /12 — hi is not an in- 
dependent parameter, the phase diagram for the general 
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FIG. 3: Finite-temperature (T = O.SSTcj ) phase diagram 
for a three-component Fermi gas with pairwise attraction be- 
tween components in the plane of chemical potential differ- 
ences /ii-/i2, Eq. (|6]). All three pairwise couplings are finite 
and such that TcjTc^ = 1.04 and Tc^/Tc^ = 0.97. As in 
the case of only two nonvanishing couplings (cf. Fig. [1} , we 
identify the regions where the normal state (N) and the su- 
perfluid states (Si) are stable. The solid segments denote 
second-order N-Si transitions. The dashed curves mark first- 
order Si-Sj transitions. Note that at the two points where 
these curves meet all three superfluids can coexist. 



case 53 ^ can be plotted in the same /ii-ft.2 plane as 
before. 



IV. DOMAIN WALL 

Until now we considered a spatially uniform system 
where one of the phases occupies the entire space. On 
the other hand, we have seen in the previous section that 
for a certain range of densities phase separation of the 
two superfluids Si and S2 occurs, as shown in Fig. 
This implies the formation of domain walls between ho- 
mogeneous phases. Similarly, domain walls must form at 
the boundaries between Si and S2 in a trapped three- 
component gas; see the text below Eq. (|25p . Let us an- 
alyze the properties of the domain wall using a grand- 
canonical thermodynamic potential. Its minima that cor- 
respond to the homogenous states Si and S2 far from the 
domain wall must be degenerate for the superfluids to 
coexist in between. Indeed, the entire phase-separated 
regions in the diagram in Fig. [2] correspond to the lines of 
degenerate minima in Fig. [T] (dashed curves). To obtain 
the domain wall solution, we need to retain the gradient 
terms in the thermodynamic potential, Eq. ([7]), and min- 



imize it subject to appropriate boundary condition. We 
first consider temperatures close to the critical one and 
later extend our considerations to lower temperatures. 



Domain walls at temperatures close to the 
critical ones 



Here, as in Sec. IIII Al we assume that conditions ([15]) 
hold. As discussed below Eq. in this case the chem- 
ical potential differences are such that \hi\/i7TT <C 1. 
Then, the prefactors in front of the gradient terms in 
Eq. ([7]) can both be approximated with 



VF 



.2^103) / 

^° 12 \2TrTJ 



V 



(32) 



We note that we cannot neglect the small differences 
of order {hi/2nT)'^ in the prefactors of the fourth-order 
terms f3i and /3i2 , as these differences enter into the equa- 
tions that determine the value of the order parameter; see 
Eq. CHI). 

For simplicity, let us assume that the translational in- 
variance is broken only along the x axis, so that the sys- 
tem is in the homogeneous state Si at x ^ —00 and 
in the state S2 at a; +00. According to Eqs. ((20|) 
and (HI]), this means 9 ^ for x —oo, 9 tt/2 
for X —t +00, and d9/dx for x — + ±00. The min- 
imization of the thermodynamic potential ([7]) yields a 
system of two second-order nonlinear differential equa- 
tion for A(a;) and 9{x) defined in Eq. (fT7)) . These equa- 
tions admit a first integral, the conserved "energy" of the 
domain wall: 



eo'(vA) 

2 



Co iV9f -I- A^ ( ai cos^ 9 + a2 sin 



1 cos^ e + (32 sin* e + 2/5i2 cos^ 9 sin^ 9 



(33) 



Our assumption p3)) implies that the two homogenous 
states Si and S2 have close values of the order parameter 
amplitude A; see Eqs. (PO)) and (PT]) . This enables us to 
neglect the (VA)^ term in Eq. ([55]) . This term changes 
little on the length scale associated with the width of the 
domain wall, while the angular variable 9{x) changes by 
tt/2 on the same length scale; i.e., the ratio of the (VA)^ 
and (V6')2 terms in Eq. §^ is of order (T^, - Tc^)/Tc,. 



Solving Eq. p9]) for A in terms of 9 and substituting the 
result into Eq. ([55]) . we arrive at 



{^oV9y 



cos 



a2 sm 



ai cos 



a2 sm 



(34) 



= -S. 



(3i cos* 9 + (32 sin* 9 + 2/?i2 cos^ 9 sin^ U 

The value of the constant S on the right-hand side can 
be determined from the boundary conditions — > and 
d9/dx — !■ as a; — !• —00: 



2(3i 



(35) 
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Note that i/S is the condensation energy density for the 
homogenous state. Indeed substituting, e.g., Eq. ([^0]) 
into Eq. (fT6|) we obtain fliy — ^ vS; see also Eq. (f3T1) . 
Using Eq. we rewrite Eq. as 



2de 
dx 



sin 29 



e VI + a cos 261 ' 



where 



and 



ai — 0(2 



(36) 



(37) 



(38) 



Here we have introduced the coherence lengths of the two 
condensates. 



and the scale factor 

Pl2 



VJW2 



1 . 



(39) 



(40) 



where a^, Pi, and /3i2 are defined by Eqs. (O, (O, and 
(llOp . respectively. In particular, to leading order in 
hi/AnT we have 



77 : 



'12 



(i) 47rr 

*(4) (1) /ll - /l; 



0.512 



47rT 

hi - h: 



(41) 



From Eq. p6p . we obtain an implicit equation for the 
spatial dependence of 6: 



X — xa 



=\/l — a, arctanh 



1 - a 



1 + a cos 29 



vT+a arctanh 



1 + a cos 26* 

1 + a 



(42) 



The parameters a and i characterize the asymmetry of 
the domain wall with respect to reflection {x — xq) 
— (x — Xq) and its width, respectively. The parameter £ 
provides a new length scale, in addition to the coherence 
lengths, via the (large) parameter ij; see Eq. ([55)1 . In the 
next subsection, we will see that the same parameter also 
enters the expression for the surface tension associated 
with the domain wall. 

An example of the spatially nonuniform order param- 
eters Ai(a;) and A2(a:) in the presence of a domain wall 
is shown in Fig. S) We also plot the angular variable 
9{x) (rescaled) and the amplitude A(x). Note that A(x) 
shows little change. This is consistent with the assump- 
tion that gradients of A(x) can be neglected near the 
critical temperature Tc^ . 

Finally, we note that because the densities on the two 
sides of the domain wall are different [see Eq. ([50]) ]. it 




FIG. 4: Profiles of the order parameter components Ai (de- 
creasing solid line) and A2 (increasing) in the presence of 
a domain wall between two superfluid states of a three- 
component Fermi gas. Here T = OmTc^, TcjTc^ = 1.05, 
and the chemical potential differences are fez/A^ = -0.66 and 
fti/A^ ~ 0.84, where T^, and A? are defined below Eq. ([TU)| . 
Note the overlap of the two components over the central re- 
gion of size I, Eq. (|38}. We also show the order parameter in 
the polar decomposition of Eq. p7|l : the dashed line is used 
for A and the dotted line for 9 x [0.063/(7r/2)]. 



could, in principle, be detected by imaging the sam- 
ple. For bosonic atoms, overlap between two Bose- 
Einstein condensates was observed long ago [2^ (for the- 
oretical studies of the two-component bosonic system, 
see Ref. [23 ). Alternatively, spatially resolved rf spec- 
troscopy [2^ could reveal the different gaps. 



B. Surface tension 

From the point of view of thermodynamic properties, 
the presence of a surface separating the two condensates 
can be taken into account by including a surface tension 
term in the thermodynamic potential |,1Q] . Moreover, as 
mentioned above, a proper treatment of surface tension 
effects is necessary to describe correctly the condensate 
profile in asymmetric traps (25j . 

The surface tension a can be calculated by integrating 
the difference between the potential in the presence of 
the domain wall (lldw) and the one in the uniform state 
(ilu) over the direction perpendicular to the domain wall: 



dx 



n, 



(43) 



Using Eq. ([5^ . we derive 



2iySe 



de 



^f{a,{P}), 



4%/^I^\/l -I- a cos 261 sin 261 
6_ cos2 26* + 2612 cos 261 + b+ 



(44) 
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with S', £, and rj defined in Eqs. 
respectively, 



M, m, and (gOD, 



(45) 



and /3i, /32, and /3i2 defined in Eqs. (O and pU]) . Using 
these definitions and hi/iirT <C 1 [see the text below 
Eq. (1221)], we estimate b+ - 0(1), ft - 0(1), 6_ - 
0{h'^/T^), and 612 - 0{h^/T^). Therefore, we replace 
the denominator in the integral (|44p by &+ and obtain 



1 r 

3a 



(46) 



By definition (|37p . the asymmetry parameter varies be- 
tween -1 and 1. Therefore, 1 > / > 23/2/3 ~ 0.94. 
Neglecting this weak dependence on the asymmetry a, 
we can write 



(47) 



This expression shows that the surface tension is deter- 
mined by the value of the condensation energy vS for the 
uniform system times the (root-mean-square) coherence 
length divided by the scale factor. As we will see in the 
next subsection, this formula for the surface tension is 
valid in a wider range of temperatures than the limiting 
case Tc, — T <C Tc, considered here. 



C. Intermediate temperatures 

In the preceding subsections we have considered a do- 
main wall near the critical temperature. On the other 
hand, as discussed at the end of Sec. lIII Al the Ginzburg- 
Landau approach remains generally valid near second- 
order phase transitions even at lower temperatures above 
Ttii; see Eq. So we can in principle analyze the prop- 
erties of the domain wall at intermediate temperatures 
(and for larger differences in the critical temperatures 
than in the previous subsections). Approaching TtH, the 
parameter /3i becomes small by definition, while in the 
superfluid state ai is finite, so we expect the difference 
between Ai and A2 to grow; see Eqs. and (HJ). If 
this is the case, the approximation in which the gradient 
of A is neglected breaks down. To remedy this, we con- 
struct in this section a variational domain wall solution. 

As a starting point for the variational approach, we 
note that the approximate domain wall solution is deter- 
mined by three parameters: the position xq, the asym- 
metry a, and the size i. The first one cannot affect the 
energy (surface tension), as it only reflects the transla- 
tional invariance of the infinite system, and henceforth 
we set Xq — 0. In the (unphysical [1^) symmetric limit 
a — > 0, we can obtain an explicit expression for, e.g., the 
profile of Ai: 



Al 



1 r 



1 — tanh 



(?) 



(48) 



This suggests the following trial functions for the order 
parameters: 



A,; 



1 T tanh ( |- T I 



(49) 



where the A^ are fixed to their asymptotic values at 
X — > ±00. We introduced a parameter 5 which describes 
the overlap between the two superfluids and enables us 
to take into account the role of the interaction term in 
Eq. ([7]). As before, we also have a parameter related 
to the domain wall thickness [3(1] . Both parameters must 
be determined by minimizing the surface tension: 



vS_ 



1 -I- coth((5) 



(50) 



with the coherence lengths 



vl Pi 
3 —ai 



which reduce to Eq. ^ as T -> Tc- . 
After minimization, cr can be written as 



(51) 



a = vS 



'7v 



with the variational scale parameter given by 
' = 5o (1 + coth(<5o)) - (1 + <5o) 



where 5ri is the solution to 



coth 5q 



So 



sinh Sq 



VWP2 
P12 



(52) 



(53) 



(54) 



Note that near the critical temperature, the right hand 
side of the above equation tends to unity, so that (5o — > 
and 77v ^ 77. Since the variational approach gives an 
upper bound on the surface tension, it also gives a lower 
one on the domain wall thickness: 



(55) 



In Fig. [5] we compare the behavior of the scale factors 
77 and ?7v as functions of temperature for T^^jT^^ = 1.04. 
The curves are computed in the two limiting cases in 
which the chemical potential differences hi are the criti- 
cal ones, hi — zkh'i{T); cf. the discussion after Eq. (|14p . 
They correspond to the points where the first-order tran- 
sition lines meet the second-order ones in Fig. [TJ There 
are two inequivalent cases depending on the relative sign 
between hi and /12. We choose these points in the phase 
diagrams because the Ginzburg-Landau expansion ([7]) is 
always valid in their vicinity as long as the temperature 
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FIG. 5: Temperature dependence of the scale factors 77V (solid 
lines) and rj (dashed lines) relating the domain wall size and 
the coherence lengths; see Eqs. (I55p and (|38|l . The coupling 
constants are chosen so that Tci/Tc2 ~ 1.04. The chemical 
potential differences hi are the critical ones; see the text after 
Eq. (|55p . They have opposite signs for the lower curves [hi — 
±hi{T), /i2 = Th2{T)] and the same sign for the upper ones 
[hi = ±hi{T), /i2 = ±/i2(r)]- In the latter case the difference 
between the two scale factors is not visible. Note that the 
scale factors are larger near and/or for chemical potential 
differences of the same sign. 

is above the tricritical temperature pS)) ; see Sec. IIIIl At 
fixed temperatures close to Tc^, we find that the scale 
parameter 77 evolves smoothly as a function of the chem- 
ical potential differences along the first-order transition 
Unes - i.e., going from one Umiting case to the other 
one. Therefore, the Hmiting cases displayed in Fig. [5] 
give upper and lower bounds on the possible values of 
the variational scale parameter ryv. 

Note that close to Tc^ the two parameters 7]^ and 77 
have very similar values, and the approximate solution 
can be trusted in this regime. At smaller temperatures 
and chemical potential differences of opposite signs, the 
approximations made in the previous sections become in- 
valid, and 77 decreases more rapidly than r]v The latter 
remains of order unity at intermediate temperatures be- 
fore quickly decreasing near T^^-^. On the contrary, for 
chemical potential differences of the same sign, both ap- 
proaches give similar results. Moreover, rj^ initially in- 
creases with decreasing temperature, leading to poten- 
tially very thick domain walls with significant overlap be- 
tween the two superfluid states. Again, when approach- 
ing Tj^j,; the scale parameter quickly decreases. How- 
ever, at these temperatures the present approach is in- 
valid - higher orders in the Ginzburg-Landau expansion 
become relevant. 

The above observations show the limits of applicability 
of the LDA (^5]) in the presence of a trap. For large rj^, 
the surface tension is small, and we expect the density 
profiles to follow the shape of the trapping potential. On 
the other hand, in this case the width of the domain 
wall is large; see Eq. (j55p . The order parameter com- 



ponents vary smoothly on this scale and density jumps 
predicted by the LDA cannot be a good approximation 
of the actual density profiles. In other words, the LDA 
breaks down, not on a length scale ^i, as usually assumed, 
but on a much longer scale. In the opposite case of small 
77v, the situation is reversed: the densities vary quickly 
on a length scale comparable to the coherence lengths. 
However, now the surface tension becomes important in 
asymmetric traps and the densities do not simply follow 
the profile of the trapping potential as in the LDA. This 
is seen, e.g., in the polarized two-component gas at low 
temperatures (25j . 

The above considerations are valid under the assump- 
tion than the sample size R is much larger than the do- 
main wall thickness £ 'sT] , in which case finite-size effects 
can be neglected. This requirement also limits the valid- 
ity of the LDA. We can estimate how large, in terms of 
the number N of trapped atoms for each species, the sam- 
ple should be in order to accommodate a domain wall. In 
the weak-coupling limit, a good estimate of the sample 
size is given by the Thomas-Fermi radius 

R ~ aU^SNy/'' , (56) 

where aho = \/ ^ I "^^ho is the harmonic oscillator length 
in the parabolic trap V{t) = ^mw^^r^. Next, we esti- 
mate the value of the coefficient ^0, Eq- ([5^ . at T ~ Tc^ 
using 

Tea 0.28£;Fe~''/2'=^l''''l , (57) 

where is the (negative) scattering length and the Fermi 
momentum (at the trap center) is 

^ —(487V) 1/6. (58) 

O-ho 

The previous three expressions ([57)l . and ([55)1 can 

be found in Ref. [ItI. Substituting Eq. ^ into Eq. (jH^ . 
we find 

^0^ 7^6^/2'=-!'^= I, (59) 
kp 

which is a lower bound for the coherence lengths defined 
in Eq. Then, for I, Eq. jHH), we can write 

^> 7^e'^/2'=^l'^=l?7. (60) 
kp 

For fcFlflsl — 1, the requirement Rji^ 1 in terms of the 
total particle number Nt ~ 3A^ becomes 

n]'"^ > 2?7 . (61) 

For a scale factor 77 ^ 10, this gives Nt ^ 10^. However, 
due to the slow growth with Nt of the left-hand side of 
Eg. (pT|) . even for a typical sample size with Nt ^ 10'' 
[28j the ratio Rjl ~ 10, and finite-size effects should be 
taken into account. Note that these estimates are also 
sensitive to the interaction strength A;/|as| [cf. Eq. ([55]) ]. 
and for a weaker interaction (e.g., A:ir|as| — 0.5) we obtain 
n\^^ > 977 instead of Eq. (HH) and Nt > 10^. 
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V. PHASE DIAGRAM AT T = 

All our previous consideration have been restricted to 
the vicinity of second-order phase transitions and hence 
to the "high-temperature" regime T > T^^.^. There is 
also a simple explicit description of the phase diagram at 
zero temperature, which we present in this section. In 
this case, all phase transitions (N-Si and S1-S2) are first 
order, and the components of the order parameter A at 
the minima of the thermodynamic potential are either 
zero or independent of the chemical potential differences 
0. Our phase diagram is in qualitative agreement with 
the numerical results of Refs. [1, [sH- Although we con- 
sider the weak-coupling regime, we expect that our re- 
sults will not qualitatively change at stronger coupling on 
the BCS side of the crossover. On the BEG side, on the 
other hand, the system behaves as a Bose-Fermi mixture 
(see, e.g., [33| for the two-component system), and we 
cannot exclude the possibility of qualitative differences 
(see also [13]). Finally, we note that in constructing the 
zero-temperature phase diagram we consider for simplic- 
ity only uniform states, neglecting the possibility that a 
spatially varying order parameter may be energetically 
favored in some regions of the phase diagram, as is the 
case for the FFLO sta te |20| in a two-component system 
dl; see, e.g., Refs. 0, Hfli, [13, HI . 

According to Eq. ([5]), the differences Sfti between the 
thermodynamic potentials in the condensed and the nor- 
mal states at the same chemical potentials is 



i = 1,2. 



(62) 



Equating 5fli to zero, we obtain the first-order transition 
lines between superfluids Si and the normal state (the 
Clogston-Chandrasekhar [3^ critical field). The condi- 
tion 6fli = dil,2 yields the first-order transition line be- 
tween the two condensates. These transitions are plotted 
as solid and dashed lines, respectively, in Fig. [6l 

Considering as before the quadratic fluctuations [cf. 
Eq. ([22]) : see also the next section], we determine the 
zero-temperature instability lines 



1,2, 



(63) 



for Si becoming unstable towards Sj. Using these ex- 
pressions, we obtain the (dotted) stability curves in the 
phase diagram shown in Fig. [S] The horizontal and verti- 
cal dotted lines indicate the instabilities of the superfluid 
states towards the normal state, which are identified by 
the conditions [^ 



2A?. 



(64) 



while the dash-dotted lines mark the instability of the 
normal state, obtained from the T ^ limit of Eq. p4)l : 



N- 



FIG. 6: Zero-temperature phase diagram for a three- 
component Fermi gas in the plane hi-h2 of chemical potential 
differences, Eq. Q. The two nonvanishing couplings con- 
stants are such that Aq/Aq = 1.05, where Aq are defined be- 
low Eq. ([5}. As at high temperature (see Fig. [l|, the normal 
state (N) is stable for large hi . Horizontal (vertical) solid seg- 
ments denote first-order N-Si (N-S2) transitions between nor- 
mal and superfiuid states (in contrast, at high temperatures 
these transitions are second order). The dashed curves iden- 
tify first-order S1-S2 transitions. The dotted curves represent 
the superfiuid-state stability limits and the dot-dashed lines 
the normal-state stability limits. Shaded areas are regions 
where both superfiuid states are (meta)stable. Note that 
these regions are much larger than the corresponding ones at 
high temperature in Fig. [T] and overlap with the normal-state 
stability regions. 



(65) 



The resulting zero-temperature phase diagram shown 
in Fig. [5] has a richer structure than that at "high" tem- 
peratures; see Fig. [TJ For example, a larger region of 
the phase diagram is occupied by metastable states due 
to the first-order nature of the transition to the normal 
state. This in turn means that in the ni-n2 density-space 
phase-separated states occupy a larger region of the phase 
diagram. Consequently, more complicated domain wall 
structures are possible that interpolate between the dif- 
ferent superfiuid states and the normal state as well, as 
is the case for a two-component Fermi gas [2^. As re- 
marked before, the phase separation translates into den- 
sity jumps in the LDA treatment of the trapping poten- 
tial. However, the validity of the LDA should be con- 
firmed by estimating the effects of domain walls and sur- 
face tension, as in the finite-temperature case. 
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VI. COLLECTIVE MODES 



Note that 



In the absence of an external potential, the existence of 
thick domain walls is a manifestation of the presence of 
soft collective modes. While the former are possible only 
in the presence of degenerate ground states, the latter 
are a more general feature of the multicomponent Fermi 
gas. In this section we present the dispersion relations 
for these modes and comment on their role in limiting the 
applicability of the BCS mean-field approach. In trapped 
Fermi gases the collective modes are known to affect the 
experimentally accessible (hydrodynamiclike) response of 
the system [40| . 

For concreteness we assume that the ground state is 
the superfluid 5*1 with homogeneous order parameter Ai 
and consider small fluctuations around this state: 

Ai (r , t) = Ai (1 + ^(r , t))e'^^'''*^ + dA2 (r , i) . (66) 

The phase fluctuations described by correspond to the 
well-known soundlike Anderson-Bogoliubov mode [4l| . 
while the amplitude fluctuations have a mass equal to 
2Ai [i^l. These two modes have also been studied in the 
BCS-BEC crossover [i^. Here we are interested in the 
fluctuations (5A2(r, t) due to pairing in the noncondensed 
channel. 

The propagator D{uj,q) of the 5A2{r,t) field is ob- 
tained by expanding Eq. ^ around the stationary point 
with Ai 7^ 0, A2 = 0: 



[i^D{L,,q)]- 



In 



A? 



(67) 



+ {u + hi - h2) H{lj, q; /ii, /i2) + J{io,q; /ii, /12) , 



[i^Do{0,0)]- 



■In 



(A?) 



/i2 — h2hi 



-lnA° = 



(70) 



yields the stability condition (|63p for i = 1. In- 
deed, since the contribution of 6A2 fluctuations to the 
action is SA2{uj,q) [D{uj,q)]~^ 6A2{uj,q), the condition 
[D{0,0)r > determines the stability of the superfluid 
Si with respect to static uniform fluctuations. 

Consider, e.g., the stability (or lack of it) of the super- 
fluid Si with respect to shifts in the chemical potentials 
in the case of equal interaction strengths. For hi = and 
small ft.2 we get 



A? 



>o 



(71) 



which shows that the superfluid Si is stable, as expected, 
since fluctuations toward condensation in the 1-3 chan- 
nel need to overcome the "Zeeman e nerg y" : cf. Eq. (|62p . 
This is contrary to the claim in Ref. [l2] that this chemi- 
cal potential shift causes the system to become unstable. 
In contrast, for h2 = the inverse propagator D{0,0)^^ 
is zero for any hi, which indicates an instability. In this 
case the stable state is the superfluid S2, as can be seen 
by repeating the above analysis with 1^2. 

Now let us determine the dispersion relation of collec- 
tive modes. For simplicity, we consider the case hi — 0. 
We have 



(72) 



where 



where functions H and J are given in the Appendix and 
A^ are the values of the zero-temperature order param- 
eter components; see the text below Eq. (jS]). Here we 
concentrate on the cases of zero temperature and vicinity 
to second-order phase transitions. Moreover, we consider 
only long wavelength fluctuations - i.e., g — > 0. 



and 



(A^)' 



hi, 



Q ( ' ' AO 



(73) 



(74) 



A. Collective modes at T = 



In the limit T ^ 0, the propagator in Eq. (I67p becomes 



[vD,,] i=ln^+7i(^) 
^2 



1 



d 



1 



hi — /i2 



n{o) 



2 2 (68) 
Vpq 



with 



H(c^) = -ln 



hi/2 



A? 



Lo + hi/2 



(69) 



(75) 



There are two branches with positive, w > 0, and neg- 
ative, u) < 0, energies. Similarly to the case of a po- 
larized normal two-component gas [4^ . we can identify 
these excitations as bifermions and biholes. We note that 
the mass of these modes explicitly depends on symmetry 
breaking due to a difference in coupling constants [flrst 
two terms on the right hand side of Eq. (|73p ] or chemical 
potentials. In the U(2)-symmetric case {hi — h2 — and 
51 = 32, so A]* = A2), the mass vanishes due to particle- 
hole symmetry. This result is independent of the weak- 
coupling assumption and holds at any coupling as long as 
particle-hole symmetry is present. Moreover, in the sym- 
metric case the collective mode speed ((7^ reduces to the 
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known result for the phase mode [4l|, vq = vp/VS. In 
other words, in the symmetric hmit in addition to the 
phase mode, there are two more modes with the same 
dispersion. This is expected in the framework of sponta- 
neous symmetry breaking from U(2) down to U(l). Due 
to condensation into the superfluid state, the system is 
invariant only under rotations that change the phase of 
the order parameter and not under rotations transform- 
ing one of the components of A = (Ai,A2,0) into the 
other. Then, to the three broken generators correspond 
three massless Goldstone bosons. On the other hand, 
in the absence of particle-hole symmetry, the dispersion 
relation is modified [l3] , and two of the massless modes 
split into a massless mode with quadratic dispersion re- 
lation and a massive one 1451 . 



B. Collective modes at finite temperatures 

Let us consider the vicinity of the second-order phase 
transition N-Si, so that Ai 0. In this case, the inverse 
propagator has a form similar to the quadratic term in 
the Ginzburg-Landau expansion ([7]) to which it reduces 
in the static limit uj ^ 0. For w 7^ the only difference is 
that the coefficients ai, (3i, and P12 depend on lo. The fre- 
quency dependence of (32 and P12 can be neglected since 
they multiply small quantities and |Aip, respectively. 



Using Eq. ([3]), we obtain 



In 

1 
2 



P2 v^q^ 



T 

1 -i{uj + h2) 

2 AttT 



/?i2|Aip 



inf 



(76) 



The general structure of this propagator is the stan- 
dard one for superconducting fluctuations [i^, with 
overdamped fluctuations typical of the time-dependent 
Ginzburg-Landau approach. What is peculiar here is 
that the mass term is proportional to |Aip. This makes 
the decay of fluctuations in the 1-3 channel (i.e., towards 
superfluid S2) faster than those toward the normal state. 
Nonetheless, they play an important role in causing devi- 
ations from mean- field theory. To show this, wc employ 
the Ginzburg-Levanyuk criterion (To| for the simple case 
/ii = /i2 = and T < • 

As is well known in the theory of second-order phase 
transitions, fluctuations strongly modify the mean-field 
behavior at temperatures close to the critical one [loj . 
The temperature window around the critical tempera- 
ture where the fluctuations dominate can be character- 
ized by the Ginzburg-Levanyuk number Gi, so that for 
£ = \T — Tc\/Tc ^ Gi fluctuations are small. In three di- 
mensions, due to amplitude fluctuations, Gi oc (Tc/Ep)^. 
This result can be obtained by writing the Ginzburg- 



Levanyuk criterion as [lO| 

TcX 



« |A|^ 



where 



X = D{0,0) oc l/ve 



is the pair susceptibility and 



f = [Dd^D-^/dq^]{0,0) oc v^p/T^ 



(77) 



(78) 



(79) 



is the coherence length squared. In both equations above 
the last term on the right is due to fluctuations of the 
order parameter Ai itself, whose propagator has the form 
similar to Eq. (|76p up to the replacement of indices 2 



1, but without the term /312IA2I ; see [46[. Using v oc 
^3/2^1/2 ^ ^ j.^^^ obtain Gi oc {TjEp)'^. 

In the present case, we can use the same approach. 
Substituting the value of the order parameter Ai, 
Eq. ((20|) . into Eq. ([76]) and using the definitions in 
Eqs. (l78|) and ([79|) . we derive for the susceptibility and 
the coherence length 



Tr- 



-1 



Xoc(i/ln^j , ^ oc UF ( TcjWlnpi 



(80) 



Using these expressions, we obtain for the three compo- 
nent case 




We see that the fluctuations in the uncondensed (1-3) 
channel shrink the region of applicability of mean-field 
theory as soon as InTc^/Tc^ ^ (Tc^/EpY - i.e., even for 
very small differences in the critical temperatures. 

In the context of the BCS-BEC crossover, we recall 
that as the strength of the interaction grows, the ra- 
tio Tc/Ep grows too. This signals the breakdown of 
the mean-field approximation as the unitary limit is 
approached from the BCS side. The above estimate 
Eq. (|8ip for the Ginzburg-Levanyuk number indicates 
that this breakdown happens much sooner in the pres- 
ence of a third interacting component. 



VII. SUMMARY AND OPEN PROBLEMS 

In this paper, we considered a three-component 
(species) Fermi gas with attractive interactions between 
fermionic species in the weak coupling regime. We con- 
firmed that there are four possible homogeneous phases: 
the normal state (N) and superfluids Si for i = 1, 2, and 
3 where species j ^ i and k ^ i,j are paired. For sim- 
plicity, for most of the paper we restricted our analysis 
to the case when the components 1 and 2 do not interact 
with each other. In this case, the homogeneous phases of 
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the system are N, Si (2 and 3 are paired), and S2 (1 and 
3 are paired). The extension of our findings to the gen- 
eral case of nonzero interaction between all components 
is straightforward; see Fig. [3] and Sec. IIII CI 

We constructed the "high" -temperature T > Ttri [see 
Eq. (fTTj)] and zero-temperature phase diagrams for arbi- 
trary differences between chemical potentials of the three 
species (Figs. [1] and HJ. In particular, we identified the 
regions where different superfluid states and the normal 
state are (meta) stable and determined the lines of first- 
order S1-S2 and second-order N-Si and N-S2 phase tran- 
sitions. We also obtained the phase diagram in the canon- 
ical ensemble in the space of particle density differences 
(713 — 77.2) and (rti — n^) (Fig. This phase diagram 
displays regions where the uniform superfiuid states are 
unstable. Phase separation between superfluids Si and 
S2 occurs for particle densities within these regions; i.e., 
the system becomes spatially inhomogeneous. 

We analyzed the properties of the domain walls be- 
tween superfluid states Si and S2. The domain walls 
are present in the phase-separated region and at an in- 
terface between layers of Si and S2 in a trapped three- 
component gas; see the text below Eq. (pSj). We de- 
termined the shape of the domain wall [see Fig. [5] and 
Eq. (|42p] and demonstrated that its thickness Eq. 
provides a new length scale that can be parametrically 
larger than the coherence lengths ^1,2 of superfluids Si, 2, 
^ ^ \/^i + • particular, this means that the two or- 
der parameters of superfluids Si and S2 can overlap sig- 
nificantly over extended regions of space. It also imposes 
severe restrictions on the LDA for evaluating the config- 
uration of superfluid and normal layers in a trap 17]; sec 
the discussion below Eq. (1^ and in the end of Sec. lIVCl 
The sharp boundaries between the superfluids predicted 
by the LDA have to be smeared over the length scale £ 
(rather than the coherence lengths ^1 or ^2 ). Further- 
more, the LDA is valid only when the size of the trap, i?, 
is much larger than the domain wall thickness, R ^ i. 
Otherwise, the two superfluids coexist throughout the 
trap. For experimentally attainable systems, the condi- 
tion R i translates into the total number of atoms 
Nt ^ 10'* with corrections to the LDA being significant 
even for typical numbers in experiments, Nt ~ 10^; see 



Eq. ([6T|) and the text after it. We also evaluated the sur- 
face tension associated with the domain wall, Eqs. ([47]) 
and (j52|) , which needs to be taken into account when con- 
sidering the shape of the interface between superfluids Si 
and S2. 

Finally, we studied the collective modes (fluctuations) 
specific to our system in Sec. I VII Namely, in the super- 
fluid state Si with order parameter Ai there are fluctu- 
ations 5/S.2{r,t) of the order parameter of superfluid S2 
and vice versa. We evaluated the mass and the dispersion 
relations of these collective modes at zero temperature 
and in the vicinity of the N-Si transition. At T = the 
mass is determined by perturbations that break the U(2) 
symmetry between species 1 and 2 - the difference in 
chemical potentials and coupling constants for the inter- 



action with 3. In the symmetric case the mass vanishes. 
At small symmetry breaking the collective modes soften 
and their mass can be parametrically smaller than the 
BCS energy gaps of superfluids Si and S2. Similarly, near 
the critical temperature of the N-Si these fluctuations 
can significantly increase the Ginzburg-Levanyuk number 
[see Eq. (j8ip ] in comparison to the two-component sys- 
tem. This indicates that stronger deviations from mean- 
field theory are possible in a three-component system. 

The results outlined above were obtained in the weak- 
coupling BCS limit. A natural question is how they are 
modified in the BCS-BEC crossover regime and in partic- 
ular at the unitary limit for two of the three components 
when the corresponding scattering length diverges. In 
the two-component case, there is a single length and en- 
ergy scale at unitarity at T = 0. This is not so in our 
case if the symmetry between the components is broken. 
Therefore, we expect qualitatively the same picture such 
as extended domain walls, soft modes, etc., as long as 
no true bound states are formed. It is also interesting to 
study these phenomena at lower temperatures close and 
below the tricritical temperature (pS)) which limits the 
applicability of our Ginzburg-Landau approach. 

Let us also emphasize that to make more quantitative 
predictions about the possible experimental realization 
and detection of coexisting multiple superfluid states, it 
is necessary to go beyond or at least improve the LDA. 
Further work is also required to understand the effects of 
fluctuations in the unpaired channel on experimentally 
accessible quantities such as critical temperatures and 
the frequencies of collective oscillations in trapped gases 
in the hydrodynamic regime [47| . 
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APPENDIX 

In this appendix we present the expressions for the 
functions H and J introduced in Eq. (p7)) for the fluctu- 
ation propagatoj\ The method to derive these functions 
is explained in 



I4I . So here we limit ourselves to the fi- 



nal results, which are straightforward extensions of those 
found in \\^ 



H 



1 



2v J (27r)3 2Ep 
1 - fiEp ~ hi 12) - /(Cp-g - hi 12 + h2) 



(82) 



uj~ Ep- £_p_g + /ii - /i2 

fiEp + hl/2)- f{Cp-g-hl/2 + h2) 



to ■ 



Ep - C 
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2v 



1 - f{Ep - - /(^p-, - + h2) 

fiEp + hi/2) - fi^p^q - h,/2 + /12) 



UJ + Ep - ^p-q + hi - h2 



where 



P 

2m 



+ M3 



(83) 



(84) 



In the hmit hi 0, Eqs. ([5^ and ([55)1 reduce (up to a 
normahzation factor) to the functions H and J obtained 



in [13. 

We note that in deriving 
the spectrum near the Fermi 
hole symmetry. Namely, we 
as p = n{pF +i/vF), where 
vp the Fermi velocity, and n 
sphere. Then, the integral 
with the integral over ^ and 



(2^)2 



, e.g., Eq. ([55]) we linearize 
surface and assume particle- 
parametrize the momentum 
Pf is the Fermi momentum, 
the unit vector on the Fermi 
over momentum is replaced 
the vector n 



(85) 



where v is the density of states at the Fermi energy. Go- 
ing beyond this approximation would enable the study of 
particle-hole asymmetry effects. 
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